Pacotes necessários:
#if(!require(PSF)) install.packages("PSF")
#if(!require(timetk)) remotes::install_github("business-science/timetk")
if(!require(EflowStats)) remotes::install_github("USGS-R/EflowStats")
pacotes <- c(
"here",
"usethis",
"data.table",
"HEobs",
# "PSF",
"tidyverse",
"lubridate",
"fs",
"checkmate",
# "xts",
# "hydroGOF",
# "ModelMetrics",
# "forecast",
# "timetk",
"EflowStats",
"hydrostats",
#"NbClust",
# "cluster",
# "cluster.datasets",
"cowplot",
# "clValid",
"ggfortify",
#"clustree",
#"dendextend",
#"factoextra",
#"FactoMineR",
#"corrplot",
#"GGally",
#"ggiraphExtra",
"kableExtra",
"tidymodels",
"fasstr"
)
# Carregar os pacotes
easypackages::libraries(pacotes)
Scripts:
source(here('R', 'load-data.R'))
source(here('R', 'utils.R'))
source(here('R', 'prep-stn-data.R'))
Os dados importados de vazão devem ser regularmente espaçados no tempo. Esta adequação das séries diárias, se necessária, pode ser realizada com a função complete_dates() do pacote {lhmetools}. Assim assegura-se que todas estações possuam 1 observação por dia e sem datas faltantes.
Metadados
data_link <- "https://www.dropbox.com/s/d40adhw66uwueet/VazoesNaturaisONS_D_87UHEsDirceuAssis_2018.dat?dl=1"
qnat_meta <- extract_metadata(NA_character_, informative = TRUE)
glimpse(qnat_meta)
## Rows: 87
## Columns: 5
## $ estacao_codigo <dbl> 18, 237, 215, 119, 190, 32, 14, 247, 1, 216, 191, 149, …
## $ latitude <dbl> -19.918751, -22.636078, -27.776389, -23.811460, -6.7985…
## $ longitude <dbl> -49.92053, -48.27274, -51.18944, -46.52390, -43.92756, …
## $ nome_estacao <chr> "A. VERMELHA", "BARRA BONITA", "BARRA GRANDE", "BILLING…
## $ municipio <chr> "A. VERMELHA", "BARRA BONITA", "BARRA GRANDE", "BILLING…
Dados
qnat_data <- qnat_dly_ons() %>%
select(date, qnat, code_stn) %>%
lhmetools::complete_dates(group = "code_stn")
glimpse(qnat_data)
## Rows: 2,796,267
## Columns: 3
## $ date <date> 1931-01-02, 1931-01-03, 1931-01-04, 1931-01-05, 1931-01-06, …
## $ code_stn <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ qnat <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
# Incluindo nome da estacao
qnat_data <- qnat_data %>%
full_join(select(qnat_meta, estacao_codigo, nome_estacao),
by = c(code_stn = "estacao_codigo")
) %>%
rename("name_stn" = nome_estacao)
Preparação dos dados para aplicação da função calc_magnifSeven com a opção de ano hidrológico (water year) para o posto 74 da ONS (G. B. Munhoz).
qnat_posto <- qnat_data %>%
sel_station(.,station = 74)
glimpse(qnat_posto)
## Rows: 32,141
## Columns: 4
## $ date <date> 1931-01-02, 1931-01-03, 1931-01-04, 1931-01-05, 1931-01-06, …
## $ code_stn <dbl> 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 7…
## $ qnat <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ name_stn <chr> "G.B. MUNHOZ", "G.B. MUNHOZ", "G.B. MUNHOZ", "G.B. MUNHOZ", "…
summary(qnat_posto)
## date code_stn qnat name_stn
## Min. :1931-01-02 Min. :74 Min. : 69.58 Length:32141
## 1st Qu.:1953-01-01 1st Qu.:74 1st Qu.: 300.02 Class :character
## Median :1975-01-01 Median :74 Median : 524.20 Mode :character
## Mean :1975-01-01 Mean :74 Mean : 730.61
## 3rd Qu.:1996-12-31 3rd Qu.:74 3rd Qu.: 958.57
## Max. :2018-12-31 Max. :74 Max. :7832.77
## NA's :13605
magnif7(stn_data = select(qnat_posto, date, qnat))
## # A tibble: 8 × 2
## indice statistic
## <chr> <dbl>
## 1 lam1 736.
## 2 tau2 0.42
## 3 tau3 0.35
## 4 tau4 0.18
## 5 ar1 0.982
## 6 amplitude 0.17
## 7 phase 262.
## 8 bfi 0.413
by_stn <- qnat_data %>%
group_by(code_stn) %>%
nest()
Testando fasstr
https://bcgov.github.io/fasstr/articles/fasstr.html
Será relevante considerar a diferenças nos períodos dos anos hidrológicos por posto?
check_season_plots <- qnat_data %>%
rename("Date" = date) %>%
plot_daily_stats(
values = "qnat",
groups = "code_stn",
start_year = 1991,
end_year = 2010,
log_discharge = TRUE,
add_year = 2001,
complete_years = TRUE,
include_title = TRUE
)
codes <- names(check_season_plots) %>%
readr::parse_number()
lookup <- qnat_data %>%
select(contains("stn")) %>%
distinct() %>%
slice(order(codes)) %>%
arrange(code_stn)
plot_l <- map(seq_along(codes),
function(i) {
# i = 1
nm <- filter(lookup, code_stn == codes[i]) %>%
pull(name_stn)
check_season_plots[[i]] + ggtitle(paste0("Posto: ", nm))
})
plot_l
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
##
## [[13]]
##
## [[14]]
##
## [[15]]
##
## [[16]]
##
## [[17]]
##
## [[18]]
##
## [[19]]
##
## [[20]]
##
## [[21]]
##
## [[22]]
##
## [[23]]
##
## [[24]]
##
## [[25]]
##
## [[26]]
##
## [[27]]
##
## [[28]]
##
## [[29]]
##
## [[30]]
##
## [[31]]
##
## [[32]]
##
## [[33]]
##
## [[34]]
##
## [[35]]
##
## [[36]]
##
## [[37]]
##
## [[38]]
##
## [[39]]
##
## [[40]]
##
## [[41]]
##
## [[42]]
##
## [[43]]
##
## [[44]]
##
## [[45]]
##
## [[46]]
## Warning: Transformation introduced infinite values in continuous y-axis
##
## [[47]]
##
## [[48]]
##
## [[49]]
##
## [[50]]
##
## [[51]]
##
## [[52]]
##
## [[53]]
##
## [[54]]
##
## [[55]]
##
## [[56]]
##
## [[57]]
##
## [[58]]
##
## [[59]]
##
## [[60]]
##
## [[61]]
##
## [[62]]
##
## [[63]]
##
## [[64]]
##
## [[65]]
##
## [[66]]
##
## [[67]]
##
## [[68]]
##
## [[69]]
##
## [[70]]
##
## [[71]]
##
## [[72]]
##
## [[73]]
##
## [[74]]
##
## [[75]]
##
## [[76]]
##
## [[77]]
##
## [[78]]
## Warning: Transformation introduced infinite values in continuous y-axis
##
## [[79]]
##
## [[80]]
##
## [[81]]
##
## [[82]]
##
## [[83]]
##
## [[84]]
##
## [[85]]
##
## [[86]]
##
## [[87]]
seven_stats <- by_stn %>%
ungroup() %>%
mutate(stats = map(data, ~.x %>% select(date, qnat) %>%magnif7))
seven_stats_tidy <- seven_stats %>%
select(stats) %>%
unnest() %>%
pivot_wider(names_from = indice, values_from = statistic) %>%
ungroup()
## Warning: `cols` is now required when using unnest().
## Please use `cols = c(stats)`
## Warning: Values are not uniquely identified; output will contain list-cols.
## * Use `values_fn = list` to suppress this warning.
## * Use `values_fn = length` to identify where the duplicates arise
## * Use `values_fn = {summary_fun}` to summarise duplicates
seven_stats_tidy
## # A tibble: 1 × 8
## lam1 tau2 tau3 tau4 ar1 amplitude phase bfi
## <list> <list> <list> <list> <list> <list> <list> <lis>
## 1 <dbl [87]> <dbl [87]> <dbl [87]> <dbl [87]> <dbl [87]> <dbl [87]> <dbl … <dbl…
saveRDS(seven_stats_tidy, file = here("output", "seven_stats.RDS"))